try this
what’s going on?
# For general stuff:
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.5.2
## ── Attaching packages ───────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.2.1 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.3
## ✓ tidyr 1.0.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.4.0
## Warning: package 'ggplot2' was built under R version 3.5.2
## Warning: package 'tibble' was built under R version 3.5.2
## Warning: package 'tidyr' was built under R version 3.5.2
## Warning: package 'purrr' was built under R version 3.5.2
## Warning: package 'dplyr' was built under R version 3.5.2
## Warning: package 'stringr' was built under R version 3.5.2
## Warning: package 'forcats' was built under R version 3.5.2
## ── Conflicts ──────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(janitor)
library(lubridate)
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
library(here)
## here() starts at /Users/juansilva/Google Drive/WINTER 2020 BREN/ESM244 ADV DATA/my-site-test
##
## Attaching package: 'here'
## The following object is masked from 'package:lubridate':
##
## here
library(paletteer)
# For ts stuff:
library(tsibble)
## Warning: package 'tsibble' was built under R version 3.5.2
##
## Attaching package: 'tsibble'
## The following objects are masked from 'package:lubridate':
##
## interval, new_interval
## The following object is masked from 'package:dplyr':
##
## id
library(fable)
## Warning: package 'fable' was built under R version 3.5.2
## Loading required package: fabletools
## Warning: package 'fabletools' was built under R version 3.5.2
library(fabletools)
library(feasts)
## Warning: package 'feasts' was built under R version 3.5.2
library(forecast)
## Warning: package 'forecast' was built under R version 3.5.2
##
## Attaching package: 'forecast'
## The following objects are masked from 'package:fabletools':
##
## GeomForecast, StatForecast
# For spatial stuff:
library(sf)
## Warning: package 'sf' was built under R version 3.5.2
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(tmap)
## Warning: package 'tmap' was built under R version 3.5.2
library(mapview)
us_renew <- read_csv(here("data", "renewables_cons_prod.csv")) %>%
clean_names()
## Parsed with column specification:
## cols(
## MSN = col_character(),
## YYYYMM = col_double(),
## Value = col_character(),
## Column_Order = col_double(),
## Description = col_character(),
## Unit = col_character()
## )
unique(us_renew$description)
## [1] "Wood Energy Production" "Biofuels Production"
## [3] "Total Biomass Energy Production" "Total Renewable Energy Production"
## [5] "Hydroelectric Power Consumption" "Geothermal Energy Consumption"
## [7] "Solar Energy Consumption" "Wind Energy Consumption"
## [9] "Wood Energy Consumption" "Waste Energy Consumption"
## [11] "Biofuels Consumption" "Total Biomass Energy Consumption"
## [13] "Total Renewable Energy Consumption"
We’ll focus on consumption data.
renew_clean <- us_renew %>%
mutate(description = str_to_lower(description)) %>%
filter(str_detect(description, pattern = "consumption")) %>%
filter(!str_detect(description, pattern = "total"))
renew_date <- renew_clean %>%
mutate(yr_mo_day = lubridate::parse_date_time(yyyymm, "ym")) %>%
mutate(month_sep = yearmonth(yr_mo_day)) %>% #coerce to tsibble `yearmonth` format
mutate(value = as.numeric(value)) %>%
drop_na(month_sep, value)
## Warning: 490 failed to parse.
## Warning: NAs introduced by coercion
# Want to parse the year and month? We may use this later...
renew_parsed <-renew_date %>%
mutate(month = month(yr_mo_day, label = TRUE)) %>%
mutate(year = year(yr_mo_day))
renew_gg <- ggplot(data = renew_date, aes(x = month_sep,
y = value,
group = description)) +
geom_line(aes(color = description))
renew_gg
renew_gg +
scale_color_paletteer_d("ghibli::LaputaMedium")
renew_ts <- as_tsibble(renew_parsed, key = description, index = month_sep)
renew_ts %>% autoplot(value)
renew_ts %>% gg_subseries(value)
# renew_ts %>% gg_season(value)
ggplot(data = renew_parsed, aes(x = month, y = value, group = year)) +
geom_line(aes(color = year)) +
facet_wrap(~description,
ncol = 1,
scales = "free",
strip.position = "right")
hydro_ts <- renew_ts %>%
filter(description == "hydroelectric power consumption")
# Explore:
hydro_ts %>% autoplot(value)
hydro_ts %>% gg_subseries(value)
# hydro_ts %>% gg_season(value)
# OK, what if gg_season() doesn't work?
# It's just a function that uses ggplot() to do things we already know how to do in ggplot()!
ggplot(hydro_ts, aes(x = month, y = value, group = year)) +
geom_line(aes(color = year))
index_by()What if we want to calculate consumption by quarter? We’ll use index_by() to tell R which “windows” to calculate a value with in.
Quarterly:
hydro_quarterly <- hydro_ts %>%
index_by(year_qu = ~ yearquarter(.)) %>% # monthly aggregates
summarise(
avg_consumption = mean(value)
)
head(hydro_quarterly)
## # A tsibble: 6 x 2 [1Q]
## year_qu avg_consumption
## <qtr> <dbl>
## 1 1973 Q1 261.
## 2 1973 Q2 255.
## 3 1973 Q3 212.
## 4 1973 Q4 225.
## 5 1974 Q1 292.
## 6 1974 Q2 290.
First, let’s check the decomposition (STL):
# Find STL decomposition
dcmp <- hydro_ts %>%
model(STL(value ~ season(window = 5)))
# View the components
# components(dcmp)
# Visualize the decomposed components
components(dcmp) %>% autoplot() +
theme_minimal()
# Let's check out the residuals:
hist(components(dcmp)$remainder)
hydro_ts %>%
ACF(value) %>%
autoplot()
MODEL FOR FORECASTING
hydro_model <- hydro_ts %>%
model(
arima = ARIMA(value),
ets = ETS(value)
) %>%
fabletools::forecast(h = "2 years")
hydro_model %>%
autoplot(filter(hydro_ts,
year(month_sep) > 2010),
level = NULL)
# Get spatial data:
world <- read_sf(dsn = here("data","TM_WORLD_BORDERS_SIMPL-0.3-1"), layer = "TM_WORLD_BORDERS_SIMPL-0.3") %>% clean_names()
# Quick & easy option to see those polygons (also for points, lines!)
mapview(world)
# ggplot (static)
world_base <- ggplot(data = world) +
geom_sf(aes(fill = pop2005),
color = NA) +
scale_fill_paletteer_c("viridis::viridis") +
theme_minimal()
world_base
# Let's crop it:
world_base +
coord_sf(xlim = c(-20, 50), ylim = c(-40, 40), expand = FALSE)